####################################
# This is code to replicate the analyses and figures from the supplementary materials
# Produces Tables 1-14, 18-19, 21
# Produces Figures 2-5, 7-9
# (Un)Natural Disasters: Electoral Cycles in Disaster Relief
# Code developed by Alicia Cooperman
####################################
# Analysis conducted using
# R version 4.0.2 (2020-06-22) 
# on Platform: x86_64-apple-darwin17.0 (64-bit)
####################################

rm(list=ls(all=TRUE))
dev.off()
library(dplyr) # Version 0.8.5
library(lfe) # Version 2.8.5
library(stargazer) # Version 5.2.2
library(sensemakr) # Version 0.1.3
library(survival) # Version 3.1.12
library(xtable) # Version 1.8.4
load("cooperman_2021_data.Rdata") 
source("cooperman_2021_functions.R") 

# Table 1: Summary Statistics ----------------------------------------------

summary <- subset(x=df, select=c(drought_bin, SPI_6_June, belowave6mo_June, temp_monthly_C_JantoJune, PET_JantoJune, drought_mayor_2years, 
                                 drought_mayor, drought_gov, lncattle, corn_share, bean_share, mayor_win_vshare_prev, 
                                 copartisan_gov, copartisan_pres, lameduck, poptot_cont, poprural_cont, percrur, inccap_cont))
cov.lab.1 <- c("Drought Declared in a Year", "SPI Jan-June","Below Average SPI Jan-June", "Temp. Jan-June", "PET Jan-June", "Drought During Mayor's Last 2 Years", "Drought in Mayor Election Year", "Drought in State/Federal Election Year", "Cattle (log)",  "Corn Share", "Beans Share","Mayor Margin Victory Previous", "Mayor Copartisan Gov", "Mayor Copartisan Pres", "Mayor Second Term",  "Total Population", "Rural Population", "Percent Rural", "Income Per Capita")
sink("SI_tab1.tex")
stargazer(summary,float=FALSE, single.row=TRUE,font.size="footnotesize",title="Summary Statistics",covariate.labels=cov.lab.1, omit.summary.stat = c("p25", "p75"), digits=2)
sink()


# Table 2: Election Years and Drought Declaration by Below Average --------

df$int_mayorelec_belowave <- df$mayorelecyear*df$belowave6mo_June

elec1 <- felm(drought_bin ~ belowave6mo_June + mayorelecyear  + govelecyear + time + time.sq  + PT  + PET_JantoJune +  lncattle + corn_share + bean_share |ibgecode|0|state.year, data=df)
elec2 <- felm(drought_bin ~ belowave6mo_June + mayorelecyear + int_mayorelec_belowave + time + time.sq  + PT + PET_JantoJune +  lncattle + corn_share + bean_share  |ibgecode|0|state.year, data=df)
elec3 <- felm(drought_bin ~ belowave6mo_June + mayorelecyear + int_mayorelec_belowave  + govelecyear*belowave6mo_June + time + time.sq + PT  + PET_JantoJune +  lncattle + corn_share + bean_share |ibgecode|0|state.year, data=df)
elec4 <- felm(drought_bin ~ belowave6mo_June + mayorelecyear + int_mayorelec_belowave + time  + PT + PET_JantoJune +  lncattle + corn_share + bean_share  |ibgecode|0|state.year, data=df)
elec5 <- felm(drought_bin ~ belowave6mo_June + mayorelecyear + int_mayorelec_belowave  + govelecyear*belowave6mo_June + time  + PT  + PET_JantoJune +  lncattle + corn_share + bean_share |ibgecode|0|state.year, data=df)
elec6 <- felm(drought_bin ~ belowave6mo_June + int_mayorelec_belowave + time + time.sq  + PT + PET_JantoJune +  lncattle + corn_share + bean_share  |ibgecode |0|state.year, data=df)
elec7 <- felm(drought_bin ~ belowave6mo_June + mayorelecyear + int_mayorelec_belowave + time + time.sq + PT + PET_JantoJune +  lncattle + corn_share + bean_share  |ibgecode + year|0|state.year, data=df)
elec8 <- felm(drought_bin ~ belowave6mo_June + mayorelecyear + int_mayorelec_belowave  + govelecyear*belowave6mo_June + time + time.sq + PT  + PET_JantoJune +  lncattle + corn_share + bean_share |ibgecode + year|0|state.year, data=df)


tab2_labels <- c("Below Average Rainfall","Mayor Election Year",  "State/Fed Election Year", "Mayor Election * Below Ave. Rainfall","Time Trend", "Time Trend Squared", "PT Party", "Temp. and Evaporation", "Cattle (log)", "Corn (share of ag. land)", "Beans (share of ag. land)", "State/Fed Election * Below Ave. Rainfall")
cat(stargazer(elec1, elec2, elec3, elec4, elec5, elec6, elec7, elec8,
              dep.var.labels="Drought Emergency Declared",
              covariate.labels=tab2_labels,float=FALSE,omit.stat=c("ll", "f"), df = F,
              add.lines=list(c("Clustered SEs", "State-Year", "State-Year", "State-Year", "State-Year", "State-Year","State-Year","State-Year", "State-Year"), c("Municipal FE", "Yes", "Yes", "Yes", "Yes", "Yes","Yes",  "Yes", "Yes"), c("Year FE", "No", "No", "No","No", "No","No", "Yes", "Yes")))
    , file="SI_tab2.tex", sep="\n")



# Table 3: Election Years and Drought Declaration by Rainfall Cate --------

elec.spicat1 <- felm(drought_bin ~ SPI_cat + mayorelecyear  + govelecyear + time  + time.sq + PT  + PET_JantoJune +  lncattle + corn_share + bean_share |ibgecode|0|state.year, data=df)
elec.spicat2 <- felm(drought_bin ~ SPI_cat*mayorelecyear + time + time.sq  + PT + PET_JantoJune +  lncattle + corn_share + bean_share  |ibgecode|0|state.year, data=df)
elec.spicat3 <- felm(drought_bin ~ SPI_cat*mayorelecyear  + SPI_cat*govelecyear + time + time.sq + PT  + PET_JantoJune +  lncattle + corn_share + bean_share |ibgecode|0|state.year, data=df)
elec.spicat4 <- felm(drought_bin ~ SPI_cat*mayorelecyear + time  + PT + PET_JantoJune +  lncattle + corn_share + bean_share  |ibgecode|0|state.year, data=df)
elec.spicat5 <- felm(drought_bin ~ SPI_cat*mayorelecyear  + SPI_cat*govelecyear + time  + PT  + PET_JantoJune +  lncattle + corn_share + bean_share |ibgecode|0|state.year, data=df)
elec.spicat6 <- felm(drought_bin ~ SPI_cat*mayorelecyear + time  + PT + PET_JantoJune +  lncattle + corn_share + bean_share  |ibgecode + year|0|state.year, data=df)
elec.spicat7 <- felm(drought_bin ~ SPI_cat*mayorelecyear  + SPI_cat*govelecyear + time + time.sq + PT  + PET_JantoJune +  lncattle + corn_share + bean_share |ibgecode + year|0|state.year, data=df)

tab3_labels <- c("Mild Drought (-1 <= SPI < 0)", "Mild Rainfall (0 <= SPI < 1)", "Very High Rainfall (1 <= SPI)",
                         "Mayor Election Year",  "State/Fed Election Year", "Time Trend", "Time Trend Squared", 
                         "PT Party", "Temp. and Evaporation", "Cattle (log)", "Corn (share of ag. land)", "Beans (share of ag. land)",
                         "Mayor Election * Mild Drought (-1 <= SPI < 0)", "Mayor Election * Mild Rainfall (0 <= SPI < 1)", "Mayor Election * Very High Rainfall (1 <= SPI)", 
                         "State/Fed Election * Mild Drought (-1 <= SPI < 0)", "State/Fed Election * Mild Rainfall (0 <= SPI < 1)", "State/Fed Election * Very High Rainfall (1 <= SPI)")
cat(stargazer(elec.spicat1, elec.spicat2, elec.spicat3, elec.spicat4, elec.spicat5, elec.spicat6, elec.spicat7,
              dep.var.labels="Drought Emergency Declared",
              float=FALSE,omit.stat=c("ll", "f"), df = F,
              covariate.labels=tab3_labels, add.lines=list(c("Clustered SEs", "State-Year", "State-Year", "State-Year", "State-Year", "State-Year","State-Year", "State-Year"), c("Municipal FE", "Yes", "Yes", "Yes", "Yes", "Yes",  "Yes", "Yes"), c("Year FE", "No", "No", "No","No", "No", "Yes", "Yes")))
    , file="SI_tab3.tex", sep="\n")


# Table 4: Election Years and Drought Declaration with Lagged Agri --------

elec.aglag  <- felm(drought_bin ~ belowave6mo_June + mayorelecyear  + govelecyear + time + time.sq + PT  + PET_JantoJune +  lncattle_lag + corn_share_lag + bean_share_lag  |ibgecode|0|state.year, data=df)
elec.may.aglag  <- felm(drought_bin ~ mayorelecyear*belowave6mo_June + time + time.sq + PT + PET_JantoJune +  lncattle_lag + corn_share_lag + bean_share_lag  |ibgecode|0|state.year, data=df)
elec.maygov.aglag  <- felm(drought_bin ~ mayorelecyear*belowave6mo_June  + govelecyear*belowave6mo_June  + time + time.sq + PT  + PET_JantoJune + lncattle_lag + corn_share_lag + bean_share_lag  |ibgecode|0|state.year, data=df)
mayorgovlag_int_labels <- c("Below Average Rainfall","Mayor Election Year",  "State/Fed Election Year","Time","Time Squared", "PT Party", "Temp. and Evaporation", "Cattle (log, lagged)", "Corn (share of ag. land, lagged)", "Beans (share of ag. land, lagged)","Mayor Election * Below Ave. Rainfall", "State/Fed Election * Below Ave. Rainfall")
cat(stargazer(elec.aglag,elec.may.aglag, elec.maygov.aglag,  dep.var.labels="Drought Emergency Declared",covariate.labels=mayorgovlag_int_labels,float=FALSE,add.lines=list(c("Clustered SEs", "State-Year", "State-Year", "State-Year")), omit.stat=c("ll", "f"))
    , file="SI_tab4.tex", sep="\n")

# Table 5: Mayoral Term and Drought Declaration ---------------------------

lameduck.below0 <- felm(drought_bin ~ lameduck + mayorelecyear + time + time.sq + PT  + PET_JantoJune +  lncattle + corn_share + bean_share |ibgecode |0|state.year, data=df, subset =SPI_6_June<0)
lameduck.below0.int <- felm(drought_bin ~ lameduck*mayorelecyear + time + time.sq + PT  + PET_JantoJune +  lncattle + corn_share + bean_share |ibgecode |0|state.year, data=df, subset =SPI_6_June<0)
lameduck.above0 <- felm(drought_bin ~ lameduck + mayorelecyear + time + time.sq + PT  + PET_JantoJune +  lncattle + corn_share + bean_share |ibgecode|0|state.year, data=df, subset =SPI_6_June>0)
lameduck.above0.int <- felm(drought_bin ~ lameduck*mayorelecyear + time + time.sq + PT  + PET_JantoJune +  lncattle + corn_share + bean_share |ibgecode|0|state.year, data=df, subset =SPI_6_June>0)
lameduck.above1 <- felm(drought_bin ~ lameduck + mayorelecyear + time + time.sq + PT  + PET_JantoJune +  lncattle + corn_share + bean_share |ibgecode|0|state.year, data=df, subset =SPI_6_June>1)
lameduck.above1.int <- felm(drought_bin ~ lameduck*mayorelecyear + time + time.sq + PT  + PET_JantoJune +  lncattle + corn_share + bean_share |ibgecode|0|state.year, data=df, subset =SPI_6_June>1)


cat(stargazer(lameduck.below0,lameduck.below0.int, lameduck.above0, lameduck.above0.int, lameduck.above1,lameduck.above1.int,
              float=FALSE, dep.var.labels="Drought Emergency Declared",
              covariate.labels = c("Second Term Mayor", "Mayor Election Year","Time", "Time Squared", "PT Party", "Temp. and Evaporation", "Cattle (log)", "Corn (share of ag. land)", "Beans (share of ag. land)", "Mayor Election * Second Term Mayor"), 
              column.labels = c("Below Ave.", "Below Ave.", "Above Ave.","Above Ave.", "Very High", "Very High"), df = F)
    , file="SI_tab5.tex", sep="\n")


# Table 6: Mayoral Term and Drought Declaration by SPI Category -  --------

df$SPI_cat <- relevel(df$SPI_cat, ref="Very high rainfall (> 1)")
lameduck.triple.felm.veryhigh <- felm(drought_bin ~ lameduck*mayorelecyear*SPI_cat + time + time.sq + PT  + PET_JantoJune +  lncattle + corn_share + bean_share |ibgecode |0|state.year, data=df)
df$SPI_cat <- relevel(df$SPI_cat, ref="Severe drought (< -1)")
lameduck.triple.felm.verylow <- felm(drought_bin ~ lameduck*mayorelecyear*SPI_cat + time + time.sq + PT  + PET_JantoJune +  lncattle + corn_share + bean_share |ibgecode |0|state.year, data=df)

covars_lameducktriple <- c("Second Term", "Mayor Election Year", "Very high rainfall (1 <= SPI)", "Severe drought (SPI < -1)", "Mild drought (-1 <= SPI < 0)", "Mild rainfall (0 <= SPI < 1)", 
                          "Time","Time Squared",  "PT Party", "Temp. and Evaporation", "Cattle (log)", "Corn (share of ag. land)", "Beans (share of ag. land)",
                           "Second Term * Mayor Election Year", "Second Term * Very high rainfall (1 <= SPI)", "Second Term * Severe drought (SPI < -1)", "Second Term * Mild drought (-1 <= SPI < 0)", "Second Term * Mild rainfall (0 <= SPI < 1)",
                           "Mayor Elec * Very high rainfall (1 <= SPI)", "Mayor Elec * Severe drought (SPI < -1)", "Mayor Elec * Mild drought (-1 <= SPI < 0)", "Mayor Elec * Mild rainfall (0 <= SPI < 1)",
                           "Second Term * Mayor Elec * Very high rainfall (1 <= SPI)", "Second Term * Mayor Elec * Severe drought (SPI < -1)", "Second Term * Mayor Elec * Mild drought (-1 <= SPI < 0)", "Second Term * Mayor Elec * Mild rainfall (0 <= SPI < 1)")
cat(stargazer(lameduck.triple.felm.verylow,lameduck.triple.felm.veryhigh, covariate.labels = covars_lameducktriple,
              float=FALSE, dep.var.labels="Drought Emergency Declared", df = F, no.space = T)
    , file="SI_tab6.tex", sep="\n")


# Table 7: Incumbent Mayor Candidacy and Drought Declaration --------------

elec.incrun <- felm(drought_bin ~ belowave6mo_June + incrun + time + time.sq + PT  + PET_JantoJune +  lncattle + corn_share + bean_share |state|0|0, data=df, subset =lameduck==0)
elec.incrun.int <- felm(drought_bin ~ belowave6mo_June*incrun + time + time.sq+PT  + PET_JantoJune +  lncattle + corn_share + bean_share |state|0|0, data=df, subset =lameduck==0)
elec.incrun.mun <- felm(drought_bin ~ belowave6mo_June + incrun +time + time.sq+ PT  + PET_JantoJune +  lncattle + corn_share + bean_share |ibgecode|0|0, data=df, subset =lameduck==0)
elec.incrun.int.mun <- felm(drought_bin ~ belowave6mo_June*incrun +time + time.sq+ PT  + PET_JantoJune +  lncattle + corn_share + bean_share |ibgecode|0|0, data=df, subset =lameduck==0)
cat(stargazer(elec.incrun,elec.incrun.int,elec.incrun.mun,elec.incrun.int.mun,
              float=FALSE, dep.var.labels="Drought Emergency Declared",
              covariate.labels = c("Below Average Rainfall", "Incumbent Runs Re-election","Time","Time Squared", "PT Party", "Temp. and Evaporation", "Cattle (log)", "Corn (share of ag. land)", "Beans (share of ag. land)", "Incumbent Runs * Below Average"),
              add.lines=list(c("FE", "State", "State", "Municipal", "Municipal")), df = F)
    , file="SI_tab7.tex", sep="\n")


# Table 8: Incumbent Mayor Candidacy and Re-Election by Relief Timing in Mayor Election Year  --------

incrun1 <- felm(incrun~ belowave6mo_June + PET_JantoJune + mayor_win_vshare_prev + copartisan_pres + copartisan_gov + lncattle + corn_share + bean_share |state + year|0|state.year, data=df, subset=lameduck==0)
incrun2 <- felm(incrun~ drought_mayor + mayor_win_vshare_prev   + copartisan_pres + copartisan_gov   + lncattle + corn_share + bean_share |state + year|0|state.year, data=df, subset=lameduck==0)
incrun3 <- felm(incrun~ drought_mayor_2years + mayor_win_vshare_prev + copartisan_pres + copartisan_gov    + lncattle + corn_share + bean_share |state + year|0|state.year, data=df, subset=lameduck==0)
incwin1 <- felm(incwin~ belowave6mo_June + PET_JantoJune + mayor_win_vshare_prev  + copartisan_pres + copartisan_gov   + lncattle + corn_share + bean_share |state + year|0|state.year, data=df, subset=lameduck==0)
incwin2 <- felm(incwin~ drought_mayor + mayor_win_vshare_prev + copartisan_pres + copartisan_gov   + lncattle + corn_share + bean_share |state + year|0|state.year, data=df, subset=lameduck==0)
incwin3 <- felm(incwin~ drought_mayor_2years + mayor_win_vshare_prev + copartisan_pres + copartisan_gov   + lncattle + corn_share + bean_share |state + year|0|state.year, data=df, subset=lameduck==0)
demandcombinedlabels <- c("Below Average Rainfall","Temp. and Evaporation", "Drought Declaration in Mayor Election Year","Drought Declaration in Last Two Years",
                          "Mayor's Vote Share in Previous Election", "Copartisan President", "Copartisan Governor","Cattle (log)", "Corn (share of ag. land)", 
                          "Beans (share of ag. land)")
cat(stargazer(incrun1, incrun2,incrun3, incwin1, incwin2,incwin3,
              covariate.labels=demandcombinedlabels,
              float=FALSE, dep.var.labels=c("Incumbent Runs","Incumbent Wins"),df = F)
    , file="SI_tab8.tex", sep="\n")


# Table 9: Incumbent Mayor Candidacy and Re-Election by Relief Timing in Mayor Election Year with Municipal FE --------


incrun1.munyearFE <- felm(incrun~ belowave6mo_June + PET_JantoJune + mayor_win_vshare_prev + copartisan_pres + copartisan_gov + lncattle + corn_share + bean_share |ibgecode + year|0|state.year, data=df, subset=lameduck==0)
incrun2.munyearFE <- felm(incrun~ drought_mayor + mayor_win_vshare_prev   + copartisan_pres + copartisan_gov   + lncattle + corn_share + bean_share |ibgecode + year|0|state.year, data=df, subset=lameduck==0)
incrun3.munyearFE <- felm(incrun~ drought_mayor_2years + mayor_win_vshare_prev + copartisan_pres + copartisan_gov    + lncattle + corn_share + bean_share |ibgecode + year|0|state.year, data=df, subset=lameduck==0)
incwin1.munyearFE <- felm(incwin~ belowave6mo_June + PET_JantoJune + mayor_win_vshare_prev  + copartisan_pres + copartisan_gov   + lncattle + corn_share + bean_share |ibgecode + year|0|state.year, data=df, subset=lameduck==0)
incwin2.munyearFE <- felm(incwin~ drought_mayor + mayor_win_vshare_prev + copartisan_pres + copartisan_gov   + lncattle + corn_share + bean_share |ibgecode + year|0|state.year, data=df, subset=lameduck==0)
incwin3.munyearFE <- felm(incwin~ drought_mayor_2years + mayor_win_vshare_prev + copartisan_pres + copartisan_gov   + lncattle + corn_share + bean_share |ibgecode + year|0|state.year, data=df, subset=lameduck==0)
demandcombinedlabels <- c("Below Average Rainfall","Temp. and Evaporation", "Drought Declaration in Mayor Election Year","Drought Declaration in Last Two Years", "Mayor's Vote Share in Previous Election",
                          "Copartisan President", "Copartisan Governor","Cattle (log)", "Corn (share of ag. land)", "Beans (share of ag. land)")
cat(stargazer(incrun1.munyearFE, incrun2.munyearFE,incrun3.munyearFE,
              incwin1.munyearFE, incwin2.munyearFE,incwin3.munyearFE,
              covariate.labels=demandcombinedlabels,
              float=FALSE, dep.var.labels=c("Incumbent Runs","Incumbent Wins"),df = F)
    , file="SI_tab9.tex", sep="\n")



# Table 10: Incumbent Mayor Candidacy with Alternative Specifications --------

incrun4 <- felm(incrun~ SPI_6_June + SPI_6_June_sq + PET_JantoJune + mayor_win_vshare_prev + copartisan_pres + copartisan_gov + lncattle + corn_share + bean_share |state + year|0|state.year, data=df, subset=lameduck==0)
incrun5 <- felm(incrun~ belowave6mo_June + PET_JantoJune + mayor_win_vshare_prev + copartisan_pres + copartisan_gov + lncattle + corn_share + bean_share |state + year|0|state.year, data=df, subset=lameduck==0)
incrun6 <- felm(incrun~ drought_mayor + mayor_win_vshare_prev   + copartisan_pres + copartisan_gov   + lncattle + corn_share + bean_share |state + year|0|state.year, data=df, subset=lameduck==0)
incrun7 <- felm(incrun~ belowave6mo_June + PET_JantoJune+ drought_mayor + mayor_win_vshare_prev   + copartisan_pres + copartisan_gov   + lncattle + corn_share + bean_share |state + year|0|state.year, data=df, subset=lameduck==0)
incrun8 <- felm(incrun~ drought_mayor_2years + mayor_win_vshare_prev + copartisan_pres + copartisan_gov    + lncattle + corn_share + bean_share |state + year|0|state.year, data=df, subset=lameduck==0)
incrun9 <- felm(incrun~ drought_mayor*belowave6mo_June + mayor_win_vshare_prev + copartisan_pres + copartisan_gov    + lncattle + corn_share + bean_share |state + year|0|state.year, data=df, subset=lameduck==0)
tab10labels <- c("SPI", "SPI-Squared",  "Below Average Rainfall","Temp. and Evaporation", "Drought Declaration in Mayor Election Year","Drought Declaration in Last Two Years",
                          "Mayor's Vote Share in Previous Election", "Copartisan President", "Copartisan Governor","Cattle (log)", "Corn (share of ag. land)", "Beans (share of ag. land)",
                          "Below Average * Declaration Mayor Elec Year")
cat(stargazer(incrun4,incrun5, incrun6,incrun7, incrun8,incrun9, covariate.labels = tab10labels,
              float=FALSE, dep.var.labels=c("Incumbent Runs"),df = F)
    , file="SI_tab10.tex", sep="\n")

# Table 11: Incumbent Mayor Wins Re-Election with Alternative Specifications --------

incwin4 <- felm(incwin~ SPI_6_June + SPI_6_June_sq + PET_JantoJune + mayor_win_vshare_prev  + copartisan_pres + copartisan_gov   + lncattle + corn_share + bean_share |state + year|0|state.year, data=df, subset=lameduck==0)
incwin5 <- felm(incwin~ belowave6mo_June + PET_JantoJune + mayor_win_vshare_prev  + copartisan_pres + copartisan_gov   + lncattle + corn_share + bean_share |state + year|0|state.year, data=df, subset=lameduck==0)
incwin6 <- felm(incwin~ drought_mayor + mayor_win_vshare_prev + copartisan_pres + copartisan_gov   + lncattle + corn_share + bean_share |state + year|0|state.year, data=df, subset=lameduck==0)
incwin7 <- felm(incwin~ belowave6mo_June + PET_JantoJune+ drought_mayor + mayor_win_vshare_prev + copartisan_pres + copartisan_gov   + lncattle + corn_share + bean_share |state + year|0|state.year, data=df, subset=lameduck==0)
incwin8 <- felm(incwin~ drought_mayor_2years + mayor_win_vshare_prev + copartisan_pres + copartisan_gov   + lncattle + corn_share + bean_share |state + year|0|state.year, data=df, subset=lameduck==0)
incwin9 <- felm(incwin~ drought_mayor*belowave6mo_June + mayor_win_vshare_prev + copartisan_pres + copartisan_gov   + lncattle + corn_share + bean_share |state + year|0|state.year, data=df, subset=lameduck==0)

tab11labels <- c("SPI", "SPI-Squared",  "Below Average Rainfall","Temp. and Evaporation", "Drought Declaration in Mayor Election Year","Drought Declaration in Last Two Years",
                          "Mayor's Vote Share in Previous Election", "Copartisan President", "Copartisan Governor","Cattle (log)", "Corn (share of ag. land)", "Beans (share of ag. land)",
                          "Below Average * Declaration Mayor Elec Year")
cat(stargazer(incwin4,incwin5, incwin6,incwin7, incwin8,incwin9, covariate.labels = tab11labels,
              float=FALSE, dep.var.labels=c("Incumbent Wins"),df = F)
    , file="SI_tab11.tex", sep="\n")




# Figures 2 and 3: Sensitivity Analysis for Incumbent Re-Election --------

incwin.lm <- lm(incwin~ drought_mayor + mayor_win_vshare_prev + copartisan_pres + copartisan_gov   + lncattle + corn_share + bean_share + as.factor(state) + as.factor(year), data=df, subset=lameduck==0)
sensitivity <- sensemakr(incwin.lm, treatment="drought_mayor", benchmark_covariates = c("mayor_win_vshare_prev"), kd=1:3)
pdf(file = "SI_fig2.pdf", onefile=FALSE,paper="special",width=5,height=5)
plot(sensitivity)
dev.off()
pdf(file = "SI_fig3.pdf", onefile=FALSE,paper="special",width=5,height=5)
plot(sensitivity, sensitivity.of="t-value")
dev.off()

# Figures 4 and 5: Sensitivity Analysis for Declarations --------


df$mayorelecyearbelowave <- df$mayorelecyear*df$belowave6mo_June
elec.may.lm <- lm(drought_bin ~ mayorelecyear + belowave6mo_June + mayorelecyearbelowave + time + time.sq + PT + PET_JantoJune +  lncattle + corn_share + bean_share + as.factor(ibgecode), data=df)
sensitivity <- sensemakr(elec.may.lm, treatment="mayorelecyearbelowave", benchmark_covariates = "PET_JantoJune", kd=1:3)
pdf(file = "SI_fig4.pdf", onefile=FALSE,paper="special",width=5,height=5)
plot(sensitivity)
dev.off()
pdf(file = "SI_fig5.pdf", onefile=FALSE,paper="special",width=5,height=5)
plot(sensitivity, sensitivity.of="t-value")
dev.off()


# Table 12: Drought Declaration with Different Partisan Alignments --------


pt.clusterstate <-felm(drought_bin ~ mayorelecyear*PT + SPI_6_June + PET_JantoJune + time + time.sq + lncattle + corn_share + bean_share |ibgecode|0|state.year, data=df, subset=PT_variation==1)
copartpres.clusterstate <- felm(drought_bin ~ mayorelecyear*copartisan_pres + SPI_6_June + PET_JantoJune + time + time.sq +lncattle + corn_share + bean_share  |ibgecode|0|state.year, data=df, subset= copartisan_pres_variation==1)
copartgov.clusterstate <- felm(drought_bin ~ mayorelecyear*copartisan_gov + SPI_6_June + PET_JantoJune + time + time.sq +lncattle + corn_share + bean_share   |ibgecode|0|state.year, data=df, subset=copartisan_gov_variation==1)
prescoal.clusterstate <- felm(drought_bin ~ mayorelecyear*pres_coalition_mayor + SPI_6_June + PET_JantoJune + time + time.sq +lncattle + corn_share + bean_share |ibgecode|0|state.year, data=df, subset= pres_coalition_mayor_variation==1)
df2 <- subset(df,!is.na(lncattle))
df2 <- subset(df2,!is.na(corn_share))
df2 <- subset(df2,!is.na(bean_share))
alignmentlabels <- c("Mayor Election Year","PT", "Copartisan Mayor President", "Copartisan Mayor Governor", "Copartisan Mayor President Coalition", "Rainfall SPI", "Temp. and Evaporation","Time","Time Squared", "Cattle (log)", "Corn (share of ag. land)", "Beans (share of ag. land)", "Mayor Elec. * PT", "Mayor Elec. * Co. Mayor Pres.", "Mayor Elect. * Co. Mayor Gov.","Mayor Elec. * Co. Mayor Pres. Coal.")
cat(stargazer(pt.clusterstate,copartpres.clusterstate,copartgov.clusterstate,prescoal.clusterstate
              ,add.lines=list(c("Number of Municipalities", length(unique(df2$ibgecode[df2$PT_variation==1])),  length(unique(df2$ibgecode[df2$copartisan_pres_variation==1])), length(unique(df2$ibgecode[df2$copartisan_gov_variation==1])), length(unique(df2$ibgecode[df2$pres_coalition_mayor_variation==1]))))
              , dep.var.labels="Drought Emergency Declared",float=FALSE, covariate.labels=alignmentlabels
              , omit.stat=c("ll", "f", "ser"), no.space=F), file="SI_tab12.tex", sep="\n")

# Table 13: Drought Declaration with Different Partisan Alignments, No Limit --------


pt.clusterstate.nolimit <-       felm(drought_bin ~ mayorelecyear*PT + SPI_6_June + PET_JantoJune + time + time.sq +lncattle + corn_share + bean_share |ibgecode|0|state.year, data=df)
copartpres.clusterstate.nolimit <- felm(drought_bin ~ mayorelecyear*copartisan_pres + SPI_6_June + PET_JantoJune + time + time.sq +lncattle + corn_share + bean_share  |ibgecode|0|state.year, data=df)
copartgov.clusterstate.nolimit <- felm(drought_bin ~ mayorelecyear*copartisan_gov + SPI_6_June + PET_JantoJune + time + time.sq +lncattle + corn_share + bean_share   |ibgecode|0|state.year, data=df)
prescoal.clusterstate.nolimit <- felm(drought_bin ~ mayorelecyear*pres_coalition_mayor + SPI_6_June + PET_JantoJune + time + time.sq +lncattle + corn_share + bean_share |ibgecode|0|state.year, data=df)


alignmentlabels <- c("Mayor Election Year","PT", "Copartisan Mayor President", "Copartisan Mayor Governor", "Copartisan Mayor President Coalition", "Rainfall SPI", "Temp. and Evaporation","Time","Time Squared", "Cattle (log)", "Corn (share of ag. land)", "Beans (share of ag. land)", "Mayor Elec. * PT", "Mayor Elec. * Co. Mayor Pres.", "Mayor Elect. * Co. Mayor Gov.","Mayor Elec. * Co. Mayor Pres. Coal.")
cat(stargazer(pt.clusterstate.nolimit,copartpres.clusterstate.nolimit,copartgov.clusterstate.nolimit,prescoal.clusterstate.nolimit
              , dep.var.labels="Drought Emergency Declared",float=FALSE, covariate.labels=alignmentlabels
              , omit.stat=c("ll", "f", "ser"), no.space=F), file="SI_tab13.tex", sep="\n")


# Table 14: Drought Declaration with PT Alignment in Below and Above Average Rainfall --------

pt.clusterstate <-felm(drought_bin ~ mayorelecyear*PT + SPI_6_June + PET_JantoJune + time + time.sq + lncattle + corn_share + bean_share |ibgecode|0|state.year, data=df, subset=PT_variation==1)
pt.clusterstate.low <- felm(drought_bin ~ mayorelecyear*PT + SPI_6_June + PET_JantoJune +time + time.sq + lncattle + corn_share + bean_share |state|0|state.year, data=df, subset=PT_variation==1 & SPI_6_June<0)
pt.clusterstate.high <- felm(drought_bin ~ mayorelecyear*PT + SPI_6_June + PET_JantoJune +time + time.sq + lncattle + corn_share + bean_share |state|0|state.year, data=df, subset=PT_variation==1 & SPI_6_June>0)

alignmentlabels <- c("Mayor Election Year","PT", "Rainfall SPI", "Temp. and Evaporation","Time", "Time Squared", "Cattle (log)", "Corn (share of ag. land)", "Beans (share of ag. land)", "Mayor Elec. * PT")
cat(stargazer(pt.clusterstate,pt.clusterstate.low,pt.clusterstate.high
              ,column.labels = c("Pooled","Below Average Rainfall", "Above Average Rainfall")
              ,add.lines=list(c("FE","Municipal", "State", "State"))
              ,dep.var.labels="Drought Emergency Declared",float=FALSE, covariate.labels=alignmentlabels
              , omit.stat=c("ll", "f", "ser"), no.space=F), file="SI_tab14.tex", sep="\n")


# Table 18: Drought Declaration with Partisan Alignment and Below Average Rainfall --------

pt.clusterstate.int <-       felm(drought_bin ~ mayorelecyear*PT*belowave6mo_June + PET_JantoJune + time + time.sq + lncattle + corn_share + bean_share |ibgecode|0|state.year, data=df, subset=PT_variation==1)
copartpres.clusterstate.int <- felm(drought_bin ~ mayorelecyear*copartisan_pres*belowave6mo_June  + PET_JantoJune + time + time.sq +lncattle + corn_share + bean_share  |ibgecode|0|state.year, data=df, subset= copartisan_pres_variation==1)
copartgov.clusterstate.int <- felm(drought_bin ~ mayorelecyear*copartisan_gov*belowave6mo_June + PET_JantoJune + time + time.sq +lncattle + corn_share + bean_share   |ibgecode|0|state.year, data=df, subset=copartisan_gov_variation==1)
prescoal.clusterstate.int <- felm(drought_bin ~ mayorelecyear*pres_coalition_mayor*belowave6mo_June+ PET_JantoJune + time + time.sq +lncattle + corn_share + bean_share |ibgecode|0|state.year, data=df, subset= pres_coalition_mayor_variation==1)

alignmentlabels <- c("Mayor Election Year","PT", "Copartisan Mayor President", "Copartisan Mayor Governor", "Copartisan Mayor President Coalition", "Below Average Rainfall", "Temp. and Evaporation","Time","Time Squared",
                     "Cattle (log)", "Corn (share of ag. land)", "Beans (share of ag. land)",
                     "Mayor Elec. * PT", "Mayor Elec. * Co. Mayor Pres.", "Mayor Elect. * Co. Mayor Gov.","Mayor Elec. * Co. Mayor Pres. Coal.",
                     "Mayor Elec. * Below Average", 
                     "PT * Below Ave." ,"PT * Below * Mayor Elec",
                     "Co. Mayor Pres. * Below Ave.","Co. Mayor Pres. * Below * Mayor Elec",
                     "Co. Mayor Gov. * Below Ave.","Co. Mayor Gov. * Below * Mayor Elec",
                     "Co. Mayor Pres. Coal. * Below Ave.","Co. Mayor Pres. Coal. * Below * Mayor Elec")
cat(stargazer(pt.clusterstate.int,copartpres.clusterstate.int,copartgov.clusterstate.int,prescoal.clusterstate.int 
              ,add.lines=list(c("Number of Municipalities", length(unique(df$ibgecode[df$PT_variation==1])),  length(unique(df$ibgecode[df$copartisan_pres_variation==1])), length(unique(df$ibgecode[df$copartisan_gov_variation==1])), length(unique(df$ibgecode[df$pres_coalition_mayor_variation==1]))))
              , dep.var.labels="Drought Emergency Declared",float=FALSE, covariate.labels=alignmentlabels
              , omit.stat=c("ll", "f", "ser"), no.space=T), file="SI_tab18.tex", sep="\n")


# Table 19: Limited Dependent Variable Model - Election Years and Drought Declared ------------------------------------------------

df$int_mayorelec_belowave <- df$mayorelecyear*df$belowave6mo_June

elec.noint <- felm(drought_bin ~ belowave6mo_June + mayorelecyear  + govelecyear + time + time.sq  + PT  + PET_JantoJune +  lncattle + corn_share + bean_share |ibgecode|0|state.year, data=df)
elec.may <- felm(drought_bin ~ belowave6mo_June + mayorelecyear + int_mayorelec_belowave + time + time.sq + PT + PET_JantoJune +  lncattle + corn_share + bean_share  |ibgecode|0|state.year, data=df)
elec.maygov <- felm(drought_bin ~ belowave6mo_June + mayorelecyear + int_mayorelec_belowave  + govelecyear*belowave6mo_June + time + time.sq + PT  + PET_JantoJune +  lncattle + corn_share + bean_share |ibgecode|0|state.year, data=df)

elec.logit <- clogit(drought_bin ~ belowave6mo_June + mayorelecyear  + govelecyear + time + time.sq + PT + PET_JantoJune + lncattle + corn_share + bean_share  + strata(ibgecode), data=df)
elec.may.logit <- clogit(drought_bin ~ belowave6mo_June*mayorelecyear + time + time.sq + PT + PET_JantoJune + lncattle + corn_share + bean_share  + strata(ibgecode), data=df)
elec.maygov.logit <- clogit(drought_bin ~ belowave6mo_June*mayorelecyear + belowave6mo_June*govelecyear + time + time.sq + PT + PET_JantoJune + lncattle + corn_share + bean_share  + strata(ibgecode), data=df)
mayor_int_labels <- c( "Below Average Rainfall", "Mayor Election Year", "State/Fed Election Year", "Time","Time Squared", "PT Party", "Temp. and Evaporation", "Cattle (log)", "Corn (share of ag. land)", "Beans (share of ag. land)","Mayor Election * Below Ave. Rainfall", "State/Fed Election * Below Ave. Rainfall")
cat(stargazer(elec.logit, elec.may.logit, elec.maygov.logit, omit=c("ibgecode"), single.row=F, float=F,
              covariate.labels=mayor_int_labels,  dep.var.labels =  c("Drought Emergency Declared"), model.names = F), 
    file="SI_tab19.tex", sep="\n")



# Figure 7: Distribution of SPI by Year across Semi-Arid Region -----------------

pdf(file = "SI_fig7.pdf", onefile=FALSE,paper="special",width=10,height=10)
boxplot(df$SPI_6_June~df$year.num, ylab="SPI Jan-June",xlab="Year",main="")
abline(h=0,col="red", lty=2)
abline(v=2000,col="red", lty=2)
dev.off()


# Figure 8: Rainfall Deviation and Drought Emergency Declaration - 3 Month SPI -------------------------------------------------------------
smoothing.pdf <- .3
bin.size.pdf <- .15

pdf(file = "SI_fig8.pdf", onefile=FALSE,paper="special",width=15,height=20)
par(mfrow=c(4,3))
statelabels <- c("AL", "BA", "CE", "MG", "PB", "PE", "PI", "RN", "SE")
statenames <- c("Alagoas (35)", "Bahia (257)", "Ceará (134)", "Minas Gerais (40)", "Paraíba (170)","Pernambuco (118)", "Piauí (109)", "Rio Grande do Norte (140)", "Sergipe (28)")
for (S in 1:length(statelabels)){
  create.rd.plot.pdf(bin.size=bin.size.pdf,smoothing=smoothing.pdf, z=df$SPI_3_June[df$state==statelabels[S] & !is.na(df$drought_bin)],y = df$drought_bin[df$state==statelabels[S] & !is.na(df$drought_bin)], state=statelabels[S], titles=statenames[S])
}
create.rd.plot.pdf(z=df$SPI_3_June,y = df$drought_bin,  bin.size=bin.size.pdf, smoothing=smoothing.pdf, titles = "ALL STATES (1031)")
dev.off()

# Figure 9: Rainfall Deviation and Drought Emergency Declaration - 9 Month SPI -------------------------------------------------------------
smoothing.pdf <- .3
bin.size.pdf <- .15

pdf(file = "SI_fig9.pdf", onefile=FALSE,paper="special",width=15,height=20)
par(mfrow=c(4,3))
statelabels <- c("AL", "BA", "CE", "MG", "PB", "PE", "PI", "RN", "SE")
statenames <- c("Alagoas (35)", "Bahia (257)", "Ceará (134)", "Minas Gerais (40)", "Paraíba (170)","Pernambuco (118)", "Piauí (109)", "Rio Grande do Norte (140)", "Sergipe (28)")
for (S in 1:length(statelabels)){
  create.rd.plot.pdf(bin.size=bin.size.pdf,smoothing=smoothing.pdf, z=df$SPI_9_June[df$state==statelabels[S] & !is.na(df$drought_bin)],y = df$drought_bin[df$state==statelabels[S] & !is.na(df$drought_bin)], state=statelabels[S], titles=statenames[S])
}
create.rd.plot.pdf(z=df$SPI_9_June,y = df$drought_bin,  bin.size=bin.size.pdf, smoothing=smoothing.pdf, titles = "ALL STATES (1031)")
dev.off()



# Table 21: Balance Table by PT Party Mayor -----------------------------------------------------------------------

balancedata <- df %>% group_by(ibgecode) %>% mutate(coast = coast[year==2011])
balancedata <- balancedata %>% group_by(ibgecode) %>% mutate(distcapital = distcapital[year==2011])
balancedata <- balancedata %>% group_by(ibgecode) %>% mutate(latitude = latitude[year==2011])
balancedata <- balancedata %>% group_by(ibgecode) %>% mutate(longitude = longitude[year==2011])
balancedata <- balancedata %>% group_by(ibgecode) %>% mutate(muniarea = muniarea[year==2011])
balancedata$d <- balancedata$PT
covs <-c("SPI_6_June", "precip_monthly_June", "PET_JantoJune", "poptot_cont","poprural_cont", "inccap_cont", "hdi_cont", "incgini_cont", "lncattle", "corn_share", "bean_share", "coast", "distcapital", "latitude", "longitude", "muniarea")
balancedata <- balancedata[,c(covs,"d")]
balancedata <- balancedata[complete.cases(balancedata),]
balancedata <- as.data.frame(balancedata)
treatset <- subset(balancedata, d==1)
controlset <- subset(balancedata, d==0)
mean_table <- sd_table <- NULL 
for(i in 1:length(covs)){
  if(class(balancedata[,covs[i]]) %in% c("integer", "numeric")){
    mean_entry <- t(as.matrix(c(mean(treatset[,covs[i]]), mean(controlset[,covs[i]]))))
    sd_entry <- t(as.matrix(c(sd(treatset[,covs[i]]), sd(controlset[,covs[i]]))))
    rownames(mean_entry) <- rownames(sd_entry) <-covs[i]
  }
  mean_table <- rbind(mean_table, mean_entry)
  sd_table <- rbind(sd_table, sd_entry)
}
b_table <- cbind(mean_table[,1], sd_table[,1], 
                 mean_table[,2], sd_table[,2],
                 mean_table[,1] - mean_table[,2],
                 sqrt(((sd_table[,1]^2)/nrow(treatset) + ((sd_table[,2])^2)/nrow(controlset)))) #eq 3.6
colnames(b_table) <- c("T Mean", "T SD", "C Mean", "C SD", "Diff", "Diff SE") 
statsig <- ifelse(abs(b_table[,5]/b_table[,6])>qnorm(.995),"***",ifelse(abs(b_table[,5]/b_table[,6])>qnorm(.975),"**",ifelse(abs(b_table[,5]/b_table[,6])>qnorm(.95),"*", " ")))
b_table <-round(b_table,3) 
b_table <- cbind(b_table,statsig)
colnames(b_table) <- c("T Mean", "T SD", "C Mean", "C SD", "Diff", "Diff SE","Stat Sig") 
nobs <- c(nrow(treatset),"", nrow(controlset),"", "", "", "")
b_table <- rbind(b_table, nobs)
rownames(b_table)<- c(covs,"N")
rownames(b_table)<- c("SPI Jan-June", "Monthly Precip. June (mm)", "PET Jan-June",
                      "Total Population", "Rural Population", "Income Per Capita", "HDI", "Income Gini",
                      "Cattle (log)",  "Corn Share", "Beans Share", "Coastal City", "Dist. to State Capital (km)",
                      "Latitude", "Longitude", "Municipal Area (sq km)", "N")
cat(print.xtable(xtable(b_table,digits=3), floating=F, include.rownames=T), file="SI_tab21.tex", sep="\n")



